#include "Preprocessor.h"
#include "Utils.h"
#include "Config.h"
#include <iostream>
#include <algorithm>
#include <cmath>
#include <map>
#include <vector>
#include <numeric>
#include <sstream>
#include <iomanip>

Preprocessor::Preprocessor() {
    // Get configuration from singleton
    auto& config = Config::getInstance();
    
    // Initialize member variables from config
    cycleSlipMethod_ = config.preprocessingSettings.cycleSlipMethod;
    weightingMode_ = config.preprocessingSettings.weightingMode;
    multipathThreshold_ = config.preprocessingSettings.multipathThreshold;
}

Preprocessor::~Preprocessor() {
    // 清理资源
}

void Preprocessor::setObservations(const std::vector<GNSSObservation>& observations) {
    observations_ = observations;
}

void Preprocessor::detectAndFixCycleSlips() {
    // Get configuration parameters
    auto& config = Config::getInstance();
    auto& prepSettings = config.preprocessingSettings;
    
    // Get maxTimeGap from config
    double maxTimeGap = prepSettings.maxTimeGap;
    bool fixCycleSlips = prepSettings.fixCycleSlips;
    
    // Group observations by satellite
    auto satGroups = groupObservationsBySatellite();
    
    // For each satellite
    for (auto& [satId, satObs] : satGroups) {
        // Sort by time
        std::sort(satObs.begin(), satObs.end(), 
                 [](const GNSSObservation& a, const GNSSObservation& b) {
                     return a.time < b.time;
                 });
        
        // Check for cycle slips
        for (size_t i = 1; i < satObs.size(); ++i) {
            bool hasSlip = false;
            
            // Check time gap - large gaps automatically considered as cycle slips
            double timeGap = satObs[i].time - satObs[i-1].time;
            if (timeGap > maxTimeGap) {
                satObs[i].hasCycleSlip = true;
                continue;
            }
            
            // Apply detection based on selected method from config
            switch (prepSettings.cycleSlipMethod) {
                case CycleSlipMethod::GeometryFree:
                    hasSlip = detectGeometryFreeCycleSlip(satObs[i-1], satObs[i], 
                                                         prepSettings.gfThreshold);
                    break;
                    
                case CycleSlipMethod::MelbourneWubbena:
                    hasSlip = detectMWCycleSlip(satObs[i-1], satObs[i], 
                                               prepSettings.mwThreshold);
                    break;
                    
                case CycleSlipMethod::TimeDifference:
                    hasSlip = detectTimeDifferenceCycleSlip(satObs[i-1], satObs[i], 
                                                           prepSettings.tdThreshold);
                    break;
                    
                case CycleSlipMethod::IonosphereRate:
                    hasSlip = detectIonosphereRateCycleSlip(satObs[i-1], satObs[i], 
                                                           prepSettings.ionoRateThreshold);
                    break;
                    
                case CycleSlipMethod::TurboEdit:
                    hasSlip = detectTurboEditCycleSlip(satObs[i-1], satObs[i]);
                    break;
                    
                case CycleSlipMethod::All:
                    // Use multiple methods
                    hasSlip = detectTurboEditCycleSlip(satObs[i-1], satObs[i]);
                    
                    if (!hasSlip) {
                        hasSlip = detectIonosphereRateCycleSlip(satObs[i-1], satObs[i], 
                                                               prepSettings.ionoRateThreshold);
                    }
                    if (!hasSlip) {
                        hasSlip = detectTimeDifferenceCycleSlip(satObs[i-1], satObs[i], 
                                                               prepSettings.tdThreshold);
                    }
                    break;
                    
                case CycleSlipMethod::None:
                default:
                    // No cycle slip detection
                    break;
            }
            
            // Mark as cycle slip if detected
            satObs[i].hasCycleSlip = hasSlip;
            
            // Fix cycle slips if configured and detected
            if (hasSlip && fixCycleSlips) {
                fixCycleSlip(satObs, i);
            }
        }
    }
    
    // Update the observations with cycle slip information
    updateObservationsFromGroups(satGroups);
    
    // Report cycle slip statistics
    size_t slipCount = std::count_if(observations_.begin(), observations_.end(),
                                   [](const GNSSObservation& obs) {
                                       return obs.hasCycleSlip;
                                   });
    
    std::cout << "Detected " << slipCount << " cycle slips out of " 
              << observations_.size() << " observations" << std::endl;
}

void Preprocessor::qualityControl(double minElevation) {
    // If no explicit elevation angle provided, use the one from config
    if (minElevation < 0) {
        auto& config = Config::getInstance();
        minElevation = config.preprocessingSettings.minElevation;
    }
    
    // 如果还没有处理过的观测数据，直接使用原始观测数据
    if (processedObservations_.empty()) {
        processedObservations_ = observations_;
    }
    
    // 基于高度角和信噪比的质量控制
    auto it = std::remove_if(processedObservations_.begin(), processedObservations_.end(),
                           [minElevation](const GNSSObservation& obs) {
                               return obs.elevation < minElevation;
                           });
    
    processedObservations_.erase(it, processedObservations_.end());
    
    // 如果使用精密星历，优化卫星位置
    if (usePreciseOrbit_ && sp3Interpolator_ != nullptr) {
        refineSatellitePositions();
    }
    
    // 移除伪距或载波相位缺失的观测
    it = std::remove_if(processedObservations_.begin(), processedObservations_.end(), 
                      [](const GNSSObservation& obs) {
                          return obs.P1 == 0.0 || obs.P2 == 0.0 || 
                                 obs.L1 == 0.0 || obs.L2 == 0.0;
                      });
    processedObservations_.erase(it, processedObservations_.end());
    
    std::cout << "质量控制完成，剩余观测数据: " << processedObservations_.size() << " 条" << std::endl;
}

std::vector<GNSSObservation> Preprocessor::getProcessedObservations() const {
    return processedObservations_;
}

bool Preprocessor::detectGeometryFreeCycleSlip(const GNSSObservation& prev, const GNSSObservation& curr, double threshold) {
    // 计算无几何组合
    double lambda1 = Utils::LIGHT_SPEED / prev.f1;
    double lambda2 = Utils::LIGHT_SPEED / prev.f2;
    
    double gf_prev = lambda1 * prev.L1 - lambda2 * prev.L2;
    double gf_curr = lambda1 * curr.L1 - lambda2 * curr.L2;
    
    // 计算相邻历元间的差值
    double delta_gf = gf_curr - gf_prev;
    
    // 判断是否有周跳
    return std::abs(delta_gf) > threshold;
}

bool Preprocessor::detectMWCycleSlip(const GNSSObservation& prev, const GNSSObservation& curr, double threshold) {
    // Melbourne-Wübbena combination eliminates both ionospheric effects and geometric distance
    // MW = [(f1*L1 - f2*L2)/(f1-f2)] - [(f1*P1 + f2*P2)/(f1+f2)]
    
    // Wavelengths
    double lambda1 = Utils::LIGHT_SPEED / curr.f1;
    double lambda2 = Utils::LIGHT_SPEED / curr.f2;
    double lambdaWL = Utils::LIGHT_SPEED / (curr.f1 - curr.f2);  // Wide-lane wavelength
    
    // Calculate wide-lane combination
    double wideLanePrev = (prev.f1 * lambda1 * prev.L1 - prev.f2 * lambda2 * prev.L2) / (prev.f1 - prev.f2);
    double wideLaneCurr = (curr.f1 * lambda1 * curr.L1 - curr.f2 * lambda2 * curr.L2) / (curr.f1 - curr.f2);
    
    // Calculate narrow-lane combination of pseudoranges
    double narrowLanePrev = (prev.f1 * prev.P1 + prev.f2 * prev.P2) / (prev.f1 + prev.f2);
    double narrowLaneCurr = (curr.f1 * curr.P1 + curr.f2 * curr.P2) / (curr.f1 + curr.f2);
    
    // Calculate Melbourne-Wübbena combination
    double mwPrev = wideLanePrev - narrowLanePrev;
    double mwCurr = wideLaneCurr - narrowLaneCurr;
    
    // Calculate the difference between epochs
    double mwDiff = std::abs(mwCurr - mwPrev);
    
    // Convert to cycles
    double mwDiffCycles = mwDiff / lambdaWL;
    
    // Compare with threshold (typically 0.5-1.0 cycles)
    if (threshold <= 0.0) {
        threshold = 0.75;  // Default value if not specified
    }
    
    return mwDiffCycles > threshold;
}

void Preprocessor::fixCycleSlip(std::vector<GNSSObservation>& satelliteObs, size_t index) {
    if (index == 0 || index >= satelliteObs.size()) {
        return;  // Cannot fix the first observation or invalid index
    }
    
    // Get a reference to the current observation but mark it as used to avoid warnings
    GNSSObservation& curr = satelliteObs[index];
    
    // Find the closest valid observations before and after
    GNSSObservation* prevValid = nullptr;
    GNSSObservation* nextValid = nullptr;
    
    // Search backward for a valid observation without cycle slip
    for (size_t i = index; i > 0; --i) {
        if (!satelliteObs[i-1].hasCycleSlip) {
            prevValid = &satelliteObs[i-1];
            break;
        }
    }
    
    // Search forward for a valid observation
    for (size_t i = index + 1; i < satelliteObs.size(); ++i) {
        if (!satelliteObs[i].hasCycleSlip) {
            nextValid = &satelliteObs[i];
            break;
        }
    }
    
    // Cannot fix if we don't have reference points
    if (!prevValid && !nextValid) {
        return;
    }
    
    // Implement actual cycle slip fixing algorithm
    // Linear interpolation method for phase observations
    if (prevValid && nextValid) {
        // Use linear interpolation between prevValid and nextValid
        double timeDiff = nextValid->time - prevValid->time;
        double timeFraction = (curr.time - prevValid->time) / timeDiff;
        
        // Get estimated change in L1 and L2 carrier phase
        double deltaL1 = nextValid->L1 - prevValid->L1;
        double deltaL2 = nextValid->L2 - prevValid->L2;
        
        // Interpolate expected values
        double expectedL1 = prevValid->L1 + deltaL1 * timeFraction;
        double expectedL2 = prevValid->L2 + deltaL2 * timeFraction;
        
        // Calculate integer cycle differences
        double cyclesDiffL1 = std::round((expectedL1 - curr.L1));
        double cyclesDiffL2 = std::round((expectedL2 - curr.L2));
        
        // Apply correction
        curr.L1 += cyclesDiffL1;
        curr.L2 += cyclesDiffL2;
        
        // Mark as fixed
        curr.hasCycleSlip = false;
    } 
    else if (prevValid) {
        // Only have previous valid point - use its values directly
        double cyclesDiffL1 = std::round((prevValid->L1 - curr.L1));
        double cyclesDiffL2 = std::round((prevValid->L2 - curr.L2));
        
        curr.L1 += cyclesDiffL1;
        curr.L2 += cyclesDiffL2;
        curr.hasCycleSlip = false;
    }
    else if (nextValid) {
        // Only have next valid point - use its values directly
        double cyclesDiffL1 = std::round((nextValid->L1 - curr.L1));
        double cyclesDiffL2 = std::round((nextValid->L2 - curr.L2));
        
        curr.L1 += cyclesDiffL1;
        curr.L2 += cyclesDiffL2;
        curr.hasCycleSlip = false;
    }
}

void Preprocessor::refineSatellitePositions() {
    if (!sp3Interpolator_ || !usePreciseOrbit_) {
        return;
    }
    
    for (auto& obs : observations_) {
        // Calculate precise satellite positions using SP3
        std::vector<double> pos(3, 0.0);
        if (sp3Interpolator_->getSatellitePosVel(obs.getSatellite(), obs.time, pos, nullptr)) {
            updateSatelliteGeometry(obs, pos);
        }
    }
}

void Preprocessor::updateSatelliteGeometry(GNSSObservation& obs, const std::vector<double>& satPos) {
    // 实现卫星几何信息计算
    if (satPos.size() < 3) {
        return;
    }
    
    // 创建从测站到卫星的向量
    std::vector<double> stationToSat = {
        satPos[0] - stationInfo_.x,
        satPos[1] - stationInfo_.y,
        satPos[2] - stationInfo_.z
    };
    
    // 计算距离(range)
    double range = std::sqrt(
        stationToSat[0] * stationToSat[0] +
        stationToSat[1] * stationToSat[1] +
        stationToSat[2] * stationToSat[2]
    );
    
    // 归一化向量
    std::vector<double> los = {
        stationToSat[0] / range,
        stationToSat[1] / range,
        stationToSat[2] / range
    };
    
    // 更新观测值的卫星位置和视线向量
    obs.satellitePosition = satPos;
    obs.lineOfSight = los;
    
    // 计算高度角和方位角
    
    // 首先将ECEF坐标转换为局部ENU坐标系
    // 获取测站的大地坐标（纬度、经度）
    double lat = stationInfo_.lat * M_PI / 180.0; // 转为弧度
    double lon = stationInfo_.lon * M_PI / 180.0; // 转为弧度
    
    // 构造从ECEF到ENU的旋转矩阵
    // [ -sin(lon)           cos(lon)          0      ]
    // [ -sin(lat)*cos(lon)  -sin(lat)*sin(lon) cos(lat) ]
    // [ cos(lat)*cos(lon)   cos(lat)*sin(lon) sin(lat)  ]
    
    // 应用矩阵变换到LOS向量
    double e = -los[0] * sin(lon) + los[1] * cos(lon);
    double n = -los[0] * sin(lat) * cos(lon) - los[1] * sin(lat) * sin(lon) + los[2] * cos(lat);
    double u = los[0] * cos(lat) * cos(lon) + los[1] * cos(lat) * sin(lon) + los[2] * sin(lat);
    
    // 计算高度角：arcsin(u)
    obs.elevation = asin(u) * 180.0 / M_PI;  // 转为角度
    
    // 计算方位角：arctan(e/n)
    obs.azimuth = atan2(e, n) * 180.0 / M_PI; // 转为角度
    if (obs.azimuth < 0) {
        obs.azimuth += 360.0;  // 确保方位角在0-360度范围内
    }
}

// 电离层变化率周跳检测
bool Preprocessor::detectIonosphereRateCycleSlip(const GNSSObservation& prev, const GNSSObservation& curr, double threshold) {
    // 计算电离层变化率
    double lambda1 = Utils::LIGHT_SPEED / prev.f1;
    double lambda2 = Utils::LIGHT_SPEED / prev.f2;
    
    // 计算前一个历元的电离层延迟
    double iono_prev = (lambda1 * prev.L1 - lambda2 * prev.L2) / (1.0 - lambda2/lambda1);
    
    // 计算当前历元的电离层延迟
    double iono_curr = (lambda1 * curr.L1 - lambda2 * curr.L2) / (1.0 - lambda2/lambda1);
    
    // 计算时间差
    double dt = curr.time - prev.time;
    if (dt <= 0) return false;
    
    // 计算电离层变化率
    double iono_rate = std::abs((iono_curr - iono_prev) / dt);
    
    // 判断是否有周跳
    return iono_rate > threshold;
}

// 时间差分周跳检测
bool Preprocessor::detectTimeDifferenceCycleSlip(const GNSSObservation& prev, const GNSSObservation& curr, double threshold) {
    // 计算L1和L2相位时间差分
    double dL1 = curr.L1 - prev.L1;
    double dL2 = curr.L2 - prev.L2;
    
    // 计算时间差
    double dt = curr.time - prev.time;
    if (dt <= 0) return false;
    
    // 计算相位变化率
    double L1_rate = dL1 / dt;
    double L2_rate = dL2 / dt;
    
    // 计算相位变化率差
    double diff_rate = std::abs(L1_rate - L2_rate);
    
    // 判断是否有周跳
    return diff_rate > threshold;
}

// 高频周跳检测
void Preprocessor::detectHighRateCycleSlips(std::vector<GNSSObservation>& satelliteObs) {
    // This method is for detecting cycle slips in high-rate data (e.g., 1Hz or higher)
    // It uses the time difference between adjacent epochs to adapt detection thresholds
    
    if (satelliteObs.size() < 3) {
        return;  // Need at least 3 observations
    }
    
    // Sort observations by time if not already sorted
    std::sort(satelliteObs.begin(), satelliteObs.end(), 
              [](const GNSSObservation& a, const GNSSObservation& b) {
                  return a.time < b.time;
              });
    
    // Get configuration parameters
    auto& config = Config::getInstance();
    auto& prepSettings = config.preprocessingSettings;
    
    // Base thresholds from config
    double baseGfThreshold = prepSettings.gfThreshold;
    double baseMwThreshold = prepSettings.mwThreshold;
    
    // Standard time difference expected for regular rate data (e.g., 30 seconds)
    double standardTimeDiff = 30.0;
    
    // Process each observation except the first one
    for (size_t i = 1; i < satelliteObs.size(); ++i) {
        // Calculate time difference with previous epoch
        double timeDiff = satelliteObs[i].time - satelliteObs[i-1].time;
        
        // Skip if time gap is too large (implies data gap, not high-rate)
        if (timeDiff > prepSettings.maxTimeGap) {
            satelliteObs[i].hasCycleSlip = true;  // Mark as cycle slip if time gap is too large
            continue;
        }
        
        // For high-rate data, scale thresholds based on time difference
        double timeRatio = timeDiff / standardTimeDiff;
        double adjustedGfThreshold = baseGfThreshold * std::sqrt(timeRatio);
        double adjustedMwThreshold = baseMwThreshold * std::sqrt(timeRatio);
        
        // Apply detection methods with adjusted thresholds
        bool hasSlip = false;
        
        // Apply selected detection method
        switch (cycleSlipMethod_) {
            case CycleSlipMethod::GeometryFree:
                hasSlip = detectGeometryFreeCycleSlip(satelliteObs[i-1], satelliteObs[i], adjustedGfThreshold);
                break;
                
            case CycleSlipMethod::MelbourneWubbena:
                hasSlip = detectMWCycleSlip(satelliteObs[i-1], satelliteObs[i], adjustedMwThreshold);
                break;
                
            case CycleSlipMethod::TurboEdit:
                // TurboEdit already handles multiple methods internally
                hasSlip = detectTurboEditCycleSlip(satelliteObs[i-1], satelliteObs[i]);
                break;
                
            default:
                // For ALL or other methods, use the most sensitive approach
                hasSlip = detectGeometryFreeCycleSlip(satelliteObs[i-1], satelliteObs[i], adjustedGfThreshold) ||
                          detectMWCycleSlip(satelliteObs[i-1], satelliteObs[i], adjustedMwThreshold);
                break;
        }
        
        satelliteObs[i].hasCycleSlip = hasSlip;
    }
}

// 多路径检测与缓解
void Preprocessor::detectAndMitigateMultipath() {
    // Calculate multipath index for each observation
    for (auto& obs : observations_) {
        obs.multipathIndex = calculateMultipathIndex(obs);
        obs.hasMultipath = obs.multipathIndex > multipathThreshold_; // Use the underscore member variable
    }
    
    // Group observations by epoch for potential spatial correlation analysis
    auto epochGroups = groupObservationsByEpoch();
    
    // Mitigate multipath effects for each satellite
    auto satGroups = groupObservationsBySatellite();
    for (auto& [satId, satObs] : satGroups) {
        if (satObs.size() > 5) {  // Need enough points for meaningful smoothing
            mitigateMultipathBySmoothing(satObs);
        }
    }
    
    // Update observations with mitigated values
    updateObservationsFromGroups(satGroups);
}

// 计算多路径指数
double Preprocessor::calculateMultipathIndex(const GNSSObservation& obs) {
    // 多路径指数计算
    double f1_sq = obs.f1 * obs.f1;
    double f2_sq = obs.f2 * obs.f2;
    double gamma = f1_sq / f2_sq;
    
    // 载波波长
    double lambda1 = Utils::LIGHT_SPEED / obs.f1;
    double lambda2 = Utils::LIGHT_SPEED / obs.f2;
    
    // 多路径组合
    double MP1 = obs.P1 - (1.0 + 2.0 / (gamma - 1.0)) * lambda1 * obs.L1 
                       + (2.0 / (gamma - 1.0)) * lambda2 * obs.L2;
    
    double MP2 = obs.P2 - (2.0 * gamma / (gamma - 1.0)) * lambda1 * obs.L1 
                       + ((2.0 * gamma) / (gamma - 1.0) - 1.0) * lambda2 * obs.L2;
    
    // 返回多路径指数（两种组合的平均绝对值）
    return (std::abs(MP1) + std::abs(MP2)) / 2.0;
}

// 使用加权平滑缓解多路径效应
void Preprocessor::mitigateMultipathBySmoothing(std::vector<GNSSObservation>& satelliteObs) {
    // 实现多路径效应缓解算法
    if (satelliteObs.size() < 3) {
        return;  // 至少需要3个观测值才能进行有意义的平滑
    }
    
    // 使用滑动窗口方法实现
    auto& prepSettings = Config::getInstance().preprocessingSettings;
    const int windowSize = prepSettings.smoothingWindow; // 从预处理设置中获取窗口大小
    
    // 验证窗口大小
    if (windowSize < 2) {
        // 窗口太小，使用指数平滑
        const double alpha = 0.3; // 平滑因子
        
        // 应用指数平滑
        for (size_t i = 1; i < satelliteObs.size(); ++i) {
            // 平滑P1
            satelliteObs[i].P1 = alpha * satelliteObs[i].P1 + (1 - alpha) * satelliteObs[i-1].P1;
            
            // 平滑P2
            satelliteObs[i].P2 = alpha * satelliteObs[i].P2 + (1 - alpha) * satelliteObs[i-1].P2;
            
            // 更新多路径指数
            satelliteObs[i].multipathIndex = calculateMultipathIndex(satelliteObs[i]);
        }
    } else {
        // 使用滑动窗口平滑
        std::vector<GNSSObservation> smoothedObs = satelliteObs; // 创建副本以存储平滑结果
        
        // 对每个观测值应用滑动窗口平滑
        for (size_t i = 0; i < satelliteObs.size(); ++i) {
            // 确定窗口的起始和结束索引
            int windowStart = std::max(0, static_cast<int>(i) - windowSize / 2);
            int windowEnd = std::min(static_cast<int>(satelliteObs.size()) - 1, static_cast<int>(i) + windowSize / 2);
            
            // 计算窗口内的权重 - 基于高度角和多路径指数
            std::vector<double> weights;
            double weightSum = 0.0;
            
            for (int j = windowStart; j <= windowEnd; ++j) {
                // 结合高度角和距离因子
                double elevationWeight = sin(satelliteObs[j].elevation * M_PI / 180.0);
                double timeDistFactor = 1.0 / (1.0 + std::abs(satelliteObs[j].time - satelliteObs[i].time));
                double weight = elevationWeight * timeDistFactor;
                
                // 如果观测值有明显的多路径，降低其权重
                if (satelliteObs[j].hasMultipath) {
                    weight *= 0.5;
                }
                
                weights.push_back(weight);
                weightSum += weight;
            }
            
            // 归一化权重
            if (weightSum > 0) {
                for (auto& w : weights) {
                    w /= weightSum;
                }
                
                // 应用加权平滑
                double smoothedP1 = 0.0;
                double smoothedP2 = 0.0;
                
                for (int j = windowStart, k = 0; j <= windowEnd; ++j, ++k) {
                    smoothedP1 += weights[k] * satelliteObs[j].P1;
                    smoothedP2 += weights[k] * satelliteObs[j].P2;
                }
                
                // 将平滑结果存储到副本中
                smoothedObs[i].P1 = smoothedP1;
                smoothedObs[i].P2 = smoothedP2;
                smoothedObs[i].multipathIndex = calculateMultipathIndex(smoothedObs[i]);
                
                // 更新多路径标志
                if (smoothedObs[i].multipathIndex < prepSettings.multipathThreshold) {
                    smoothedObs[i].hasMultipath = false;
                }
            }
        }
        
        // 用平滑后的结果替换原始观测值
        satelliteObs = smoothedObs;
    }
}

// 为观测值计算权重
void Preprocessor::computeObservationWeights() {
    for (auto& obs : observations_) {
        // Default weight
        obs.weight = 1.0;
        
        switch (weightingMode_) {
            case WeightingMode::Elevation: {
                // Elevation-based weighting
                obs.weight = computeElevationWeight(obs.elevation);
                break;
            }
            case WeightingMode::Snr: {
                // SNR-based weighting
                obs.weight = computeSNRWeight(obs.snr);
                break;
            }
            case WeightingMode::Multipath: {
                // Multipath-based weighting
                obs.weight = computeMultipathWeight(obs.multipathIndex);
                break;
            }
            case WeightingMode::Combined: {
                // Combined weighting strategy
                double elevWeight = computeElevationWeight(obs.elevation);
                double snrWeight = computeSNRWeight(obs.snr);
                double mpWeight = computeMultipathWeight(obs.multipathIndex);
                
                // Combine weights (e.g., multiply or average)
                obs.weight = elevWeight * snrWeight * mpWeight;
                break;
            }
            case WeightingMode::Uniform:
            default:
                // Uniform weighting (already set to 1.0)
                break;
        }
        
        // Ensure minimum weight
        const double MIN_WEIGHT = 0.01;
        if (obs.weight < MIN_WEIGHT) {
            obs.weight = MIN_WEIGHT;
        }
    }
    
    std::cout << "计算了 " << observations_.size() << " 条观测记录的权重" << std::endl;
}

// 根据高度角计算权重
double Preprocessor::computeElevationWeight(double elevation) {
    // 基于sin^2(elevation)模型
    double elev_rad = elevation * Utils::DEG2RAD;
    return std::sin(elev_rad) * std::sin(elev_rad);
}

// 根据信噪比计算权重
double Preprocessor::computeSNRWeight(double snr) {
    // 简化的SNR权重模型
    if (snr <= 0) return 0.1;  // 最小权重
    
    // 线性映射：SNR从20到50映射到权重0.1到1.0
    double weight = (snr - 20.0) / 30.0;
    weight = std::max(0.1, std::min(1.0, weight));
    
    return weight;
}

// 根据多路径指数计算权重
double Preprocessor::computeMultipathWeight(double multipathIndex) {
    // 简化的多路径权重模型：多路径指数越大，权重越小
    double maxMP = 2.0;  // 最大可接受多路径指数
    
    double norm_mp = std::abs(multipathIndex) / maxMP;
    norm_mp = std::min(1.0, norm_mp);
    
    return 1.0 - 0.9 * norm_mp;  // 映射到0.1-1.0
}

// 辅助方法
std::map<std::string, std::vector<GNSSObservation>> Preprocessor::groupObservationsByEpoch() {
    std::map<std::string, std::vector<GNSSObservation>> epochGroups;
    
    for (const auto& obs : observations_) {
        // Convert time to string with sufficient precision to identify epochs
        std::ostringstream timeStr;
        timeStr << std::fixed << std::setprecision(3) << obs.time;
        epochGroups[timeStr.str()].push_back(obs);
    }
    
    return epochGroups;
}

std::map<int, std::vector<GNSSObservation>> Preprocessor::groupObservationsBySatellite() {
    std::map<int, std::vector<GNSSObservation>> satGroups;
    
    for (const auto& obs : observations_) {
        satGroups[obs.satellite].push_back(obs);
    }
    
    // Sort each group by time
    for (auto& [satId, satObs] : satGroups) {
        std::sort(satObs.begin(), satObs.end(), 
                 [](const GNSSObservation& a, const GNSSObservation& b) {
                     return a.time < b.time;
                 });
    }
    
    return satGroups;
}

bool Preprocessor::detectTurboEditCycleSlip(const GNSSObservation& prev, const GNSSObservation& curr) {
    // Get configuration values
    auto& config = Config::getInstance();
    auto& prepSettings = config.preprocessingSettings;
    
    // Get thresholds from config
    double gfThreshold = prepSettings.gfThreshold;
    double mwThreshold = prepSettings.mwThreshold;
    
    // Detect using both methods
    bool gfDetected = detectGeometryFreeCycleSlip(prev, curr, gfThreshold);
    bool mwDetected = detectMWCycleSlip(prev, curr, mwThreshold);
    
    // TurboEdit strategy has 3 options:
    // 1. CONSERVATIVE: Detect if either method detects (OR logic)
    // 2. STRICT: Detect only if both methods detect (AND logic)
    // 3. ADAPTIVE: Use either AND or OR based on data quality
    
    // Default to conservative approach
    bool isSlip = gfDetected || mwDetected;
    
    // If we have a strict policy, require both methods to detect
    if (prepSettings.turboEditStrategy == Config::PreprocessingSettings::TurboEditStrategy::Strict) {
        isSlip = gfDetected && mwDetected;
    }
    // If using adaptive approach, choose based on data quality
    else if (prepSettings.turboEditStrategy == Config::PreprocessingSettings::TurboEditStrategy::Adaptive) {
        double qualityMetric = calculateDataQuality(prev, curr);
        if (qualityMetric >= prepSettings.adaptiveQualityThreshold) {
            isSlip = gfDetected && mwDetected;  // Strict for good quality data
        } else {
            isSlip = gfDetected || mwDetected;  // Conservative for poor quality data
        }
    }
    
    return isSlip;
}

// Helper method to calculate data quality for adaptive strategy
double Preprocessor::calculateDataQuality(const GNSSObservation& prev, const GNSSObservation& curr) {
    // Mark prev as used to avoid warning
    (void) prev;
    
    // Simple quality metric based on signal strength and elevation
    double elevationFactor = curr.elevation / 90.0; // Normalize to 0-1
    elevationFactor = std::max(0.1, elevationFactor); // Min 0.1
    
    double snrFactor = std::min(1.0, curr.snr / 50.0); // Normalize SNR to 0-1
    
    // Combined quality metric (0.0-1.0)
    return 0.6 * elevationFactor + 0.4 * snrFactor;
}

void Preprocessor::updateObservationsFromGroups(const std::map<int, std::vector<GNSSObservation>>& groups) {
    // Clear observations and rebuild from updated groups
    observations_.clear();
    
    // Append all observations from each group
    for (const auto& [key, groupObs] : groups) {
        observations_.insert(observations_.end(), groupObs.begin(), groupObs.end());
    }
    
    // Sort observations by time
    std::sort(observations_.begin(), observations_.end(), 
              [](const GNSSObservation& a, const GNSSObservation& b) {
                  return a.time < b.time;
              });
    
    std::cout << "Updated observations from groups, total: " << observations_.size() << std::endl;
}

void Preprocessor::setCycleSlipMethod(CycleSlipMethod method) {
    cycleSlipMethod_ = method;
}

void Preprocessor::setWeightingMode(WeightingMode mode) {
    weightingMode_ = mode;
}

void Preprocessor::setMultipathThreshold(double threshold) {
    multipathThreshold_ = threshold;
} 